%% 用遗传算法来优化pid参数
clear
clc
close all
%% 基本设置
file = 'danji';
load_system([file,'.slx']);%加载simulink模块
lb = [0.9 0.12 0.5];%pid参数变动的下界
ub = [1.1 0.08 0.7];%pid参数变动的上界

npop=100;%种群个数
nvar=3;%参数量
maxit=50;%迭代次数
pc1=0.8;%交叉率
pc2=0.6;
pm1=0.2;%变异率
pm2=0.1;

rng(0);%固定随机数

template.x=[];
template.y=[];
template.dominationset=[];%支配集，指个体能够支配的其他个体所对应的下标
template.dominated=[];%被支配数，指个体在自然选择中被支配的次数
template.rank=[];%等级，指个体在自然选择中所占据的生态位
template.cd=[];%拥挤度
 
pop=repmat(template,npop,1);
%% 交叉变异迭代
lhs_rand = lhsdesign(npop,nvar);
for i=1:npop
    tic
    pop(i).x = lb+lhs_rand(i,:).*(ub-lb);%随机生成参数组合
    pop(i).y=respond(file,pop(i).x);%根据参数组合得到对应的y值
    toc
end
figure('Name','遗传算法迭代过程')%迭代过程
y_best = [];
for it=1:maxit
    tic
    npc=1;
    popc=repmat(template,npop/2,2);
    fall = [pop.y];%fall存储所有的y值
    fall = reshape(fall,[2, npop]);%第一行存储所有的y1，第二行存储所有的y2
    fmax =2; %这一代中综合y值最高的个体，以其作为归一化的上限
    fave =mean(fall(1,:))./max(fall(1,:))+mean(fall(2,:))./max(fall(2,:)); %这一代中综合y值的平均数
    for i=1:npop/2%交叉
        ind=randperm(npop,2);%选择交叉个体

        fcross1=pop(ind(1)).y(1)./max(fall(1,:))+pop(ind(1)).y(2)./max(fall(2,:));
        fcross2=pop(ind(2)).y(1)./max(fall(1,:))+pop(ind(2)).y(2)./max(fall(2,:));
        fcross = max([fcross1 fcross2]);
        pc=pc2;
        if(fcross>fave)
            pc = pc1-(pc1-pc2).*(fcross-fave)./(fmax-fave);
        end

        value = rand();
        if(value<=pc)
        [popc(npc,1).x,popc(npc,2).x]=Cross(pop(ind(1)).x,pop(ind(2)).x);%交叉得到新的参数组合
        popc(npc,1).y=respond(file,popc(npc,1).x);  
        popc(npc,2).y=respond(file,popc(npc,2).x);
        npc=npc+1;
        end
    end

    npc = npc-1;
    popc(npc+1:npop/2,:)=[];
    npm=1;
    popm=repmat(template,npop,1);
    i_mutate = [];%存储变异个体的编号
    for j=1:npop%变异
        ind=randperm(npop,1);%选择变异个体

        fmutate=pop(ind(1)).y(1)./max(fall(1,:))+pop(ind(1)).y(2)./max(fall(2,:));
        pm=pm2;
        if(fmutate>fave)
            pm = pm1-(pm1-pm2).*(fmutate-fave)./(fmax-fave);
        end
        value=rand();
        if(value<=pm)
        i_mutate = [i_mutate ind];
        end
    end
    mutate_rand = lhsdesign(length(i_mutate),1);
    for j = i_mutate
        popm(npm,1).x=Mutate(pop(ind).x,lb,ub);
        popm(npm,1).y=respond(file,popm(npm).x);
        npm=npm+1;
    end

    npm=npm-1;
    popm(npm+1:npop)=[];
    popc=popc(:);

    newpop=[pop;popc;popm];%新种群
    [newpop,F]=Non_dominate_sort(newpop);
    newpop=Crowd(newpop,F);%计算拥挤度
    newpop=nsga2Sort(newpop);%根据生态等级与拥挤度进行排序
    pop=newpop(1:npop);%只取前npop个个体存活至下一次迭代

    y1=zeros(1,npop);
    y2=zeros(1,npop);
    ys=[pop.y];
    for j=1:npop
        y1(j)=ys(2*j-1);
        y2(j)=ys(2*j);
    end
    plot(1./y1,1./y2,'r*');%y1是超调量，y2是调节时间
    numtitle=num2str(it);
    title('迭代次数=',numtitle);
    xlabel('y1');
    ylabel('y2');
    set(gcf,'color','white');
    pause(0.001);
    %frame=getframe(gcf);
    %writeVideo(v,frame);%记录迭代过程
    toc
    y_sum = 1*fall(1,:)+1*fall(2,:);%
    y_best(it) = max(y_sum);
    i_best = find(y_sum == y_best(it));
end